/*
* DenseSvd.cs
* 
* Copyright (c) 2005, dnAnalytics. All rights reserved.
*/

#region Using Directives
using System;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.Math;
using dnAnalytics.Exceptions;
using dnAnalytics.Resources;
#endregion

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    /// <summary>
    /// The Compute method locks the object and call InternalCompute. InternalCompute will 
    /// decompose the matrix, set the values of u,v,s, and rank. It will will also set the 
    /// initial matrix to null.
    /// </summary>
    internal sealed class DenseSvd : Algorithm, ISvd
    {
        #region Fields
        private readonly bool computeVectors;
        private bool converged;
        private readonly int rows;
        private readonly int columns;
        private DenseMatrix matrix;
        private DenseMatrix u;
        private DenseMatrix v;
        private DenseVector s;
        private int rank;
        #endregion Fields

        #region Constructors
        /// <summary>
        /// Constructs an QR object for the given matrix.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="computeVectors">Whether to compute the left and right singular vectors or not.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <b>null</b>.</exception>
        public DenseSvd(Matrix matrix, bool computeVectors)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix", Strings.NullParameterException);
            }

            this.matrix = new DenseMatrix(matrix);
            this.computeVectors = computeVectors;
            this.rows = matrix.Rows;
            this.columns = matrix.Columns;
        }

        #endregion Constructors

        #region Public Methods

        /// <summary>
        /// Returns whether the decomposition converged or not.
        /// </summary>
        /// <returns><b>true</b> if the decomposition converged; <b>false</b> otherwise.</returns>
        public bool Converged()
        {
            Compute();
            return converged;
        }

        ///<summary>Returns the two norm of the <see cref="Matrix"/>.</summary>
        ///<returns>The two norm of the <see cref="Matrix"/>.</returns>        
        public double Norm2()
        {
            return s[0];
        }

        ///<summary>Returns the condition number <b>max(S) / min(S)</b>.</summary>
        ///<returns>The condition number.</returns>
        public double ConditionNumber()
        {
            int tmp = System.Math.Min(rows, columns) - 1;
            return s[0] / s[tmp];
        }

        ///<summary>Returns the effective numerical matrix rank.</summary>
        ///<returns>The number of non-negligible singular values.</returns>
        public int Rank()
        {
            return rank;
        }

        ///<summary>Returns the left singular vectors as a <see cref="Matrix"/>.</summary>
        ///<returns>The left singular vectors. The matrix will be <b>null</b>,
        ///if <b>computeVectors</b> in the constructor is set to false.</returns>
        public Matrix U()
        {
            if (computeVectors)
            {
                return u.Clone();
            }
            else
            {
                return null;
            }
        }

        /// <summary>
        /// Copies the left singular vectors into the <paramref name="result"/> <see cref="Matrix"/>.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the left singular vectors into.</param>
        /// <remarks>The method does nothing if <b>computeVectors</b> in the constructor is set to false.</remarks>
        public void U(Matrix result)
        {
            if (computeVectors)
            {
                u.CopyTo(result);
            }
        }

        ///<summary>Returns the right singular vectors as a <see cref="Matrix"/>.</summary>
        ///<returns>The right singular vectors. The matrix will be <b>null</b>,
        ///if <b>computeVectors</b> in the constructor is set to false.</returns>
        ///<remarks>This is the transpose of the V matrix.</remarks>
        public Matrix VT()
        {
            if (computeVectors)
            {
                return v.Clone();
            }
            else
            {
                return null;
            }
        }

        /// <summary>
        /// Copies the right singular vectors into the <paramref name="result"/> <see cref="Matrix"/>.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the right singular vectors into.</param>
        /// <remarks>The method does nothing if <b>computeVectors</b> in the constructor is set to false.</remarks>
        ///<remarks>This is the transpose of the V matrix.</remarks>
        public void VT(Matrix result)
        {
            if (computeVectors)
            {
                v.CopyTo(result);
            }
        }

        ///<summary>Returns the singular values as a diagonal <see cref="Matrix"/>.</summary>
        ///<returns>The singular values as a diagonal <see cref="Matrix"/>.</returns>        
        public Matrix W()
        {
            Matrix result = MatrixBuilder.CreateMatrix(rows, columns);
            W(result);
            return result;
        }

        /// <summary>
        /// Resets the <paramref name="result"/> <see cref="Matrix"/> and copies the singular values onto the diagonal.
        /// </summary>
        /// <param name="result"><see cref="Matrix"/> to copy the singular values to.</param>
        public void W(Matrix result)
        {
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < columns; j++)
                {
                    if (i == j)
                    {
                        result.ValueAt(i, i, s[i]);
                    }
                }
            }
        }

        ///<summary>Returns the singular values as a <see cref="Vector"/>.</summary>
        ///<returns>the singular values as a <see cref="Vector"/>.</returns>
        public Vector S()
        {
            return s.Clone();
        }

        /// <summary>
        /// Copies the singular values to the given <see cref="Vector"/>.
        /// </summary>
        /// <param name="result">The <see cref="Vector"/> to copy the singular values into.</param>
        public void S(Vector result)
        {
            for (int i = 0; i < s.Count; i++)
            {
                result[i] = s[i];
            }
        }

        /// <summary>
        /// Finds the least squares solution of <b>A*X = B</b>, where <b>A.Rows &gt;= A.Columns</b>
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <returns>The left hand side <see cref="Matrix"/>, <b>X</b>.</returns>
        /// <exception cref="InvalidOperationException">The left and right singular values were not computed, 
        /// but they are required for this operation.</exception>
        public Matrix Solve(Matrix rightHandSide)
        {
            if (u == null)
            {
                throw new InvalidMatrixOperationException(Strings.SingularVectorsNotComputed);
            }

            Matrix result = MatrixBuilder.CreateMatrix(columns, rightHandSide.Columns, MatrixType.Dense);
            Solve(rightHandSide, result);
            return result;
        }

        /// <summary>
        /// Finds the least squares solution of <b>A*X = B</b>, where <b>A.Rows &gt;= A.Columns</b>
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Matrix"/>, <b>B</b>.</param>
        /// <param name="result">The left hand side <see cref="Matrix"/>, <b>X</b>.</param>
        /// <exception cref="InvalidOperationException">The left and right singular values were not computed, 
        /// but they are required for this operation.</exception>
        public void Solve(Matrix rightHandSide, Matrix result)
        {
            if (u == null)
            {
                throw new InvalidMatrixOperationException(Strings.SingularVectorsNotComputed);
            }

            int mn = System.Math.Min(rows, columns);
            int bn = rightHandSide.Columns;

            Vector tmp = VectorBuilder.CreateVector(rows, VectorType.Dense);
            double value;

            for (int k = 0; k < bn; k++)
            {
                for (int j = 0; j < rows; j++)
                {
                    value = 0;
                    if (j < mn)
                    {
                        for (int i = 0; i < rows; i++)
                        {
                            value += u[i, j] * rightHandSide[i, k];
                        }
                        value /= s[j];
                    }
                    tmp[j] = value;
                }

                for (int j = 0; j < columns; j++)
                {
                    value = 0;
                    for (int i = 0; i < columns; i++)
                    {
                        value += v[i, j] * tmp[i];
                    }
                    result[j, k] = value;
                }
            }
        }

        /// <summary>
        /// Finds the least squares solution of <b>A*x = b</b>, where <b>A.Rows &gt;= A.Columns</b>
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Vector"/>, <b>b</b>.</param>
        /// <returns>The left hand side <see cref="Vector"/>, <b>x</b>.</returns>
        /// <exception cref="InvalidOperationException">The left and right singular values were not computed, 
        /// but they are required for this operation.</exception>
        public Vector Solve(Vector rightHandSide)
        {
            if (u == null)
            {
                throw new InvalidMatrixOperationException(Strings.SingularVectorsNotComputed);
            }

            Vector result = VectorBuilder.CreateVector(columns, VectorType.Dense);
            Solve(rightHandSide, result);
            return result;
        }

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>.
        /// </summary>
        /// <param name="rightHandSide">The right hand side <see cref="Vector"/>, <b>b</b>.</param>
        /// <param name="result">The left hand side <see cref="Vector"/>, <b>x</b>.</param>
        /// <exception cref="InvalidOperationException">The left and right singular values were not computed, 
        /// but they are required for this operation.</exception>
        public void Solve(Vector rightHandSide, Vector result)
        {
            if (u == null)
            {
                throw new InvalidMatrixOperationException(Strings.SingularVectorsNotComputed);
            }

            int mn = System.Math.Min(rows, columns);

            Vector tmp = VectorBuilder.CreateVector(rows, VectorType.Dense);
            double value;
            for (int j = 0; j < rows; j++)
            {
                value = 0;
                if (j < mn)
                {
                    for (int i = 0; i < rows; i++)
                    {
                        value += u[i, j] * rightHandSide[i];
                    }
                    value /= s[j];
                }
                tmp[j] = value;
            }

            for (int j = 0; j < columns; j++)
            {
                value = 0;
                for (int i = 0; i < columns; i++)
                {
                    value += v[i, j] * tmp[i];
                }
                result[j] = value;
            }
        }

        #endregion Public Methods

        #region Protected Members
        /// <summary>
        /// Delegates the actual computation to the correct QR subtype.
        /// </summary>
        protected override void InternalCompute()
        {
            int nm = System.Math.Min(rows + 1, columns);
            s = (DenseVector)VectorBuilder.CreateVector(nm, VectorType.Dense);

            int info = 0;
            if (computeVectors)
            {
                u = (DenseMatrix)MatrixBuilder.CreateMatrix(rows, MatrixType.Dense);
                v = (DenseMatrix)MatrixBuilder.CreateMatrix(columns, MatrixType.Dense);
                info = matrix.lapack.SVDDecomposition(rows, columns, matrix.data, s.data, u.data, v.data);
            }
            else
            {
                info = matrix.lapack.SVDDecomposition(rows, columns, matrix.data, s.data, null, null);
            }

            if (info == 0)
            {
                converged = true;
            }

            //adjust the size of s if row < columns.
            //we are using ported copy of linpack's svd code and it uses
            //a singular vector of length rows+1 when rows < columns.
            //the last element is not used and needs to be removed.
            //we should port lapack's svd routine to remove this problem.
            if (rows < columns)
            {
                nm--;
                DenseVector tmp = (DenseVector)VectorBuilder.CreateVector(nm, VectorType.Dense);
                for (int i = 0; i < nm; i++)
                {
                    tmp[i] = s[i];
                }
                s = tmp;
            }

            double eps = System.Math.Pow(2.0, -52.0);
            double tol = System.Math.Max(rows, columns) * s[0] * eps;
            rank = 0;
            for (int h = 0; h < nm; h++)
            {
                if (s[h] > tol)
                {
                    rank++;
                }
            }

            //we no longer need the original matrix.
            matrix = null;
        }
        #endregion Protected Members
    }
}
